Simulation of Soluble and Bound VEGF-stimulated in vitro Capillary-like Network Formation on Deformed Substrate

Capillary plexus cultivation is crucial in tissue engineering and regenerative medicine. Theoretical simulations have been conducted to supplement the expensive experimental works. However, the mechanisms connecting mechanical and chemical stimuli remained undefined, and the functions of the different VEGF forms in the culture environment were still unclear. In this paper, we developed a hybrid model for simulating short-term in vitro capillary incubations. We used the Cellular Potts model to predict individual cell migration, morphology change, and continuum mechanics to quantify biogel deformation and VEGF transport dynamics. By bridging the mechanical regulation and chemical stimulation in the model, the results showed good agreement between the predicted network topology and experiments, in which elongated cells connected, forming the network cords and round cells gathered, creating cobblestone-like aggregates. The results revealed that the capillary-like networks could develop in high integrity only when the mechanical and chemical couplings worked adequately, with the cell morphology and haptotaxis driven by the soluble and bound forms of VEGF, respectively, functioning simultaneously.


Introduction
Capillary development, such as in angiogenesis and vasculogenesis, plays a crucial role in embryonic development and wound healing, as well as the realization of regenerative medicine and tissue engineering.Cultivating blood capillaries is vital for the circulation connection and waste removal between artificial tissue and the living body [1][2][3].Capillary cultivation in vitro can be classified into two categories depending on the experimental settings and lasting times [4].The short-term cultivations typically last from a few hours to days.The necessary chemicals during the capillary developing process are offered well before the trial startup.Such tests usually aim to research the endothelial cell migration behavior or the interactions between the cells and their environments [5][6][7].The long-term incubations, on the other hand, are generally over weeks or months to allow for the extracellular matrix (ECM) or the functional proteins secreted from cells, mimicking the realistic environments for capillary growth and showing potential for applications such as wound healing and transplanting [8][9][10].
Since the 1980s, researchers have tried various substrates in in vitro capillary incubating experiments.They studied the relationships between endothelial cell behaviors and the mechanical properties of biogel materials.These studies revealed biogel remodeling, in which the cells exert traction forces on the substrates and disturb the material density during cell migration [5,9,11].For different purposes, researchers have designed various biogels composed of specific components, such as the mixing of collagen and gelatin to simulate the environments where proteins are denatured in burns and scalds [12] or the coated fibronectin on scaffolds to study the fibronectin structures and functions for cell attachment promotion [13,14].Fibronectin is a biomaterial widely studied in recent years for its integrin-binding sub-domains, which can bind the integrin on the cellular membrane to stabilize the cell attachment.Moreover, together with its growth factor-binding sub-domains corresponding to various growth factors, fibronectin regulates growth factor sources in both in vitro and in vivo experiments [15][16][17][18].
Later on, thanks to the development of protein purification techniques, the mainstream research topics shifted to the endothelial cell responses induced by vascular endothelial growth factors (VEGFs).The VEGF research can be traced back to 1977 when the vascular permeability factor (VPF) was discovered around the tumor cells to increase the capillary permeability and obtain extra nutrients [19].Studies have confirmed that VEGFs present multiple functions for regulating endothelial cell growth, proliferation, migration, and morphology.VEGFs comprise five isoforms: the A-, B-, C-, and D-type and the placental growth factor (PGF), which induce the differentiation of endothelial cells into blood or lymphatic vessels [20,21].Among the five isoforms, the VEGF-A, comprising the most abundant 165-amino acids in the human body, is most related to capillary development [22].Furthermore, the VEGF binding to substrates is commonly called the bound VEGF, which is considered to dominate the angiogenesis and vasculogenesis [23][24][25].Besides, experimental results in recent years showed that the synergy between integrin and bound VEGF on fibronectin enhances the VEGF stimulation on endothelial cells, further promoting angiogenesis [26].
In addition to experiments, theoretical models and numerical simulations have been conducted to study capillary development.The simulation models have two categories: continuum-and cell-based models, according to their calculation methods for cell behaviors.The continuum-based model uses partial differential equations (PDEs) to describe the cell behaviors in a continuous spatial distribution.Murray and Oster first established an equation set that includes the dynamics of cell colony and extracellular matrix (ECM).Considering the force balance between the cells and ECM, they performed the linear stability analysis to explore the cellular network conditions [27].Manoussaki et al. followed and simplified Murray and Oster's equations into two-dimensional by considering the migration of endothelial cells limited on the upper layer of the matrix, successfully obtaining the mechanical force-induced visualized pattern of non-uniform distributed cells [28].Namy et al. added the long-ranged elasticity descriptions for the matrix fibers to the theoretical model, making the calculation results more quantitatively comparable between experiments and numerical simulations [29].In addition, Gamba et al., in the era of growing understanding of VEGF, assumed VEGF to be the driver of cell migration and obtained the chemical signal-induced cellular networks [30].
Cell-based models, such as the cellular Potts model (CPM), were developed to consider cell individual behaviors using a probability approach rather than that of the cell colony.The CPM framework, viewing individual cells moving on a fixed grid, can be traced back to Graner and Glazier, who extended the Potts model with several constraints suitable for living cells, and the cell alliances tended to be stable [31,32].Then, Merks et al. coupled both the cell-based and the continuum-based models.In their hybrid approach, CPM predicted the endothelial cell migration, while PDEs calculated the VEGF concentration [33].Scianna and Munaron further included the intracellular calcium signaling induced by VEGF in their model, linking the VEGF to the cytoskeletal remodeling and elongation phenomena of endothelial cells [34].Later, the VEGF sources were modified from the autocrine assumption to paracrine based on experimental evidence.Ko ¨hn-Luque et al. assumed that haptotaxis induced by the bound form of VEGF dominated the cell migration rather than the chemotaxis generated by the soluble form of VEGF [35].Lima et al. combined the macro tumor with the micro sprouting models, suggesting that VEGF para-secreted by tumor cells provided endothelial proliferation [36].
On the other hand, van Oers et al. published a mechanical-based hybrid model with their successful simulations of cell networks and sprouts, showing that the capillary-like structures could develop via the pure stress and strain mechanisms between cells and substrates [37].After considering several phenomena, such as vascular adaptation, vessel remodeling, and capillary collapsing, Vilanova et al. developed a model based on cell movement, expansion, and apoptosis, simulating the capillary growth consistent with the experimental data [38].In recent years, Phillips et al. combined the confocal experimental data with computational simulations, predicting the total sprouting length via calibrating and verifying endothelial cell growth cycles [39,40].
The complete driving mechanisms underlying the capillary development remain unclear.There was a consensus that the bound VEGF located on the substrates dominates the developing; however, the role of the soluble VEGF in the culture medium was commonly ignored [23][24][25]35].Most mathematical models included chemical or mechanical stimulation in the development separately; the integrated connection between the mechanical and chemical mechanisms was still absent [29,41,42].Besides, although the computational results met experiments in certain conditions, several problems remain.For instance, mechanical mechanisms omitted the necessity of bound VEGF [42], whereas chemical mechanisms could not fully reproduce the results of the tilted biogel experiments [43].
This paper studies how mechanical regulation and chemical stimulation simultaneously work in capillary network development.We developed a hybrid model using cell-based CPM to predict the endothelial cell migration and morphology and using continuum-based PDEs to quantify the biogel deformation and VEGF transport dynamics.The simulation settings and the capillary-like network results were referred mainly to the short-term experiments in the literature [44] without considering endothelial cells secreting and intaking VEGF.Close attention was paid to the impact of VEGFs in the solute and bound form, respectively.In the following, we first simulate an experimental case to ensure the correction of the computation in Section 2. We show how endothelial cell elongation induced by soluble VEGF and haptotaxis driven by bound VEGF regulate the network-forming process.We then performed parameter analysis for the network formation under various mechanical and chemical mechanisms.Discussions on the model features and limitations are in Section 3, and the summarized conclusions are in Section 4. Finally, we specify the mathematical formulation and computational procedure in Section 5.

Cellular network formation
Previous experimental studies showed capillary-like networks form at 18 hours of incubation [44].Fig 1 compares the networks of the reference result in the literature and our present simulation.Table 1 presents the simulation parameters for their typical values and citations.The overall network topology agreed well.The local cell alliances were similar, too, exhibiting two regional cell colony features.First, the elongated cells connected head-to-tail and lined up, creating the network cords.Second, the remaining flattened, round cells gathered closely, forming a cobblestone-like aggregate.[7,45,46].
Fig 2 shows the computed cell distribution, the biogel density, and the bound VEGF concentration through the network formation process.Because of the diffusion effect, the simulated soluble VEGF was virtually uniform across the Petri dish.We thus omit to show the soluble VEGF results.The cell network evolution and features at different times were consistent with the experimental data [44].Initially, the cells spread randomly, and the biogel and soluble VEGF were uniform in space.The cells migrated and formed into the cellular cord elements after the first 6 hours.Then, the cellular cord elements were connected, and the cords enclosed the cell-free lacunae through 6 to 24 hours of incubation.The cell patterns altered less after 18 hours, which indicates a stable network had developed.The overall cellular network was fully developed with 1 to 3 cells in cord width, as shown in Figs 1B and 2A, which agreed with the experimental observations and could meet the sequential cord hollowing requirement in angiogenesis and vasculogenesis [44,52].
The results present a strong correlation between the cells, the biogel, and the VEGF fields, which formed a positive feedback loop among the cells, biogel, and bound VEGF.Once the endothelial cells were attached to the biogel surface, the gel, as pulled by the cell traction, started concentrating, as observed during the first few Monte-Carlo steps in the simulation.In the next stage, the concentrated biogel led to the concentrated bound VEGF; the higher the biogel density, the more the VEGF binding sits.Due to haptotaxis driven by bound VEGF, the cells moved toward the more bound VEGF areas, leading to further concentrated biogel.As the process continued, the elevation of biogel density and bound VEGF concentration would finally reach their quasi-equilibrium levels owing to the material deformation limitations and the force equilibrium, respectively, and the quasi-steady co-localization of cells, biogel, and bound VEGF appeared.

Effects of soluble and bound VEGF
To investigate the roles of the two VEGF forms in the network developing process, we compare different network features by numerically turning on/off the cellular morphological change and haptotaxis effects.For comparison, the numbers of network junctions and segments and the total segment length from the reference experiment [44]  would show a significantly lower value in junction number than in the segment number and total segment length because the coincident contact between the elongated cells would generate segments.
We also analyzed the simulated cell shapes.As shown in Fig 5 , we first evaluated the roundness index and aspect ratio for 400 cells every six incubation hours under morphology-on and haptotaxis-on.The results showed the increasing features of the elongated cells.The number of cells with a roundness index less than 0.3 grew significantly over time, as well as an aspect ratio greater than 3. Fig 6 shows the roundness and aspect ratio data collected at 24 hours for various morphology/haptotaxis turned on/off.The cells remained round if the morphology effect was absent; almost no cell presented a roundness index of less than 0.3, and more than 75% of cells had an aspect ratio between 1 and 2. No significant difference in cell shape occurred between the haptotaxis on and off cases, indicating that the morphology change cue from soluble VEGF was the primary cause for cell shape distribution.

Modality of cell migration
The discrete CPM could trace the movement of each cell.The cell volume center was recorded in each MCS to indicate its moving direction, and its orientation through time was grouped into 12 direction intervals for statistical analysis (Fig 7A).We recorded the accumulated steps a cell moved in certain direction intervals.The dip test was conducted to determine the moving direction modality for a cell during numerical cultivation.The uniform and the unimodal distributions were tested sequentially using α 1 = α 2 = 0.1 [53].Fig 7B shows the three typical migration modalities, where the isotropic mode refers to the cell performing pure random walks with no preferred directions, the one-direction mode represents the cell likely to migrate toward a specific direction significantly in addition to random walks, and the head-and-tail means the cell migrating mainly along the two opposite directions.The results agreed with the experiment observation for the elongated cells [7,54,55].We counted the number of cells in different migration modalities on a two-hour base.The classification among 400 cells in each time slot is shown in Fig 8 .For the two cases having the haptotaxis effect, the number of cells moving in one-direction mode was slightly more significant during the first few hours due to the cell aggregation from separated to the cord structures.For the two cases allowing cell morphology change, the cells appeared in the head-andtail trend at 6 hours until the cells got elongated.It was worth mentioning that about 350 cells were showing no directional favors because of the intense random walks in the first few hours and the stable network development in the end.Approximately 40 cells appeared in the onedirection mode, and the cells in the head-and-tail mode were the fewest.

VEGF parameter analysis
Five parameters related to the VEGFs in the simulation are the cell morphology change coefficient μ m , the haptotaxis coefficient μ h , the initial solute concentration c 0 , the VEGF binding rate k on and unbinding rate k off .The first two parameters regulate the strength of the chemical cues, and the other three balance the two forms of VEGF contents.The total number of junctions was used to present the cellular network complexity.The junction number is better than the segment number and length in assessing the network integrity, as shown in Fig 4 , because the elongated cells might isolated in the cavities, which should not contribute to the network integrity; however, they would be counted in the segment number and length.Fig 9 shows how the junction number at 24 hours changed for varying one of the five parameters.We conducted regression fits based on the observed data trends.The junction number at 24 hours approached an asymptotic value with increasing parameters.Thus, the regression curves take  similar intercepts, slopes, and plateau values (Fig 9D and 9E).The regression curve f(t) proposed based on the observed data trends reads the following form  The fitted values of coefficients are listed in Table 3 with the coefficient of determination r 2 = 0.9472 and adjusted coefficient of determination r 2 = 0.9447.The first exponential term denotes the junction number of the final stage, and the regressions suggested that increasing the haptotaxis strength μ h and the soluble VEGF amount c 0 raised the junction number more than expanding the morphology strength μ m under the same magnification rates.The term in the square bracket accounts for the time rate of the network development.The rising speed of the junction number was higher if we multiplied the initial soluble VEGF concentration c 0 .Indeed, increasing the initial VEGF amount was equivalent to increasing both the soluble and bound VEGF concentration contemporarily, so it played a more significant effect than other parameters.

Biogel parameter analysis
As for the parameters related to biogel deformation in the simulation, we focused on interactive parameters for the three parts: cells, biogel, and VEGFs.The first was the strength of cell traction force κ, which regulates the biogel responses from cell adhesion.The second was the ratio of VEGF binding sites γ, which bridges the communications between the VEGF and biogel.Finally, the biogel thickness h 0 was also considered for verification purposes.Other parameters related to biogel properties were not changed in value in investigations.
Fig 11 shows the junction number at 24 hours.We changed these parameters around their typical values and performed the regression fits based on the data trends.The regression equation relating the junction number with the strength of cell traction force κ is The regression equation for the junction number with the ratio of VEGF binding sites γ is The last set of data points for the various biogel thicknesses yields the logistic function as follows The fitted values of coefficients are listed in Table 4.The results showed positive correlations between the network junction and the three parameters.Eq (3) reveals the greater traction force resulted in higher biogel deformation.Eq (4) shows that increasing the VEGF binding sites led to increased bound VEGF accumulation for inducing cell migration.Eq (5) presents that the predicted junction number approached an asymptotic value as the gel increased in thickness, for the upper surface of the thicker biogel received less restriction from the Petri dish.All these results depicted the correlation between the biogel deformation and  the cellular network formation, which agreed with the experimental observations [43].In summary, the capillary-like network development relied on adequate biogel deformation, a sufficient amount of bound VEGF, and the gradient of bound VEGF to offer spatial cues for the endothelial cells.The cellular network could be well-developed only when the mechanical and chemical coupling worked simultaneously.

Discussion
Our simulations showed that the cellular-like network was well-developed if the mechanical and chemical coupling mechanisms could work simultaneously.We compare the cellular networks between the literature experiment [44] and the current simulation at 18 hours (Fig 1), even though we have collected the data until 24 hours.The experiment showed no significant drop in the number of living cells before the last 6 hours, supporting the apoptosis-neglected assumption for simulation, which is valid till 18 hours.Literature has shown that after 18 hours of incubation, the hollowed spaces inside the cellular cords could be observed [7,43].
We thus regarded the stable capillary-like network at 18 hours as our simulation target.The agreement between simulation and experiment results presents the reliability of our mathematical model.Our results successfully simulated the formation of cellular networks and the co-existence of cobblestone alliances observed in previous experiments [7,44,56].Based on experimental observation, we suppressed the cell elongation if it was surrounded by other cells in the simulation.Detailed mechanisms behind this phenomenon require further investigation.
The morphologized endothelial cells from previous studies showed the following features: elongation in shape, cytoskeletal solidification, and cytoplasm merging between connected cells [46].The arching of the elongated cells benefited the subsequent cord hollowing, so a slight decrease in the occupied top-view area was observed under the phase-contrast microscope [46,52].Although our modeling only involved the elongation probability for each cell, the solidification of the cell cytoskeleton was partially presented even in the current 2-D simulation.Our cellular shape constraints for area and perimeter were balanced in the CPM.Thus, smaller cells were obtained compared to the original.However, Our simulation ignored the cytoplasm merging between the cells.This stabilization issue was negligible in normal conditions when cells gathered due to haptotaxis but would lead to mismatches if the haptotaxis was absent.Repeated attachment and detachment of elongated cells were observed in the simulation case of morphology-on and haptotaxis-off (Fig 3C).In contrast, an incomplete but maintainable network in the experiments without haptotaxis driven by bound VEGF was shown since the cytoplasm merge stabilized the connected cells [44].
The present model focused on the cell elongation cued by the soluble VEGF without considering the detailed intracellular events activated by the soluble VEGF.The chemical-based model of Scianna and Munaron (2011) considered a series of intracellular events, including activating protein kinase and releasing arachidonic acid, nitric oxide, and calcium ions [34].Their model may provide a basis for further investigating the intracellular signal path for cell elongation cued by soluble VEGF.The VEGF secretion from and later consumption by endothelial cells was also included in their model, which resulted in a non-constant total VEGF amount.The VEGF consumption could decrease the VEGF concentration over time, reducing cellular morphology and haptotaxis tendencies, and might cause an earlier but less complete network [42].
Our model considered the bound VEGF was accumulated along with the biogel around the cells through the cell traction, which is consistent with the finding of Malandrino et al. ( 2019) [57].They showed that cells can accumulate ECM around them through pulling.In our simulation, the concentrated bound VEGF under the cells drove the haptotaxis and the alignment of the elongated cells.However, the present model has ignored the transport of soluble VEGF hindered by the attached cells on the substrate.The VEGF diffusion and binding and unbinding rates on the cell-attached surface were as if on a cell-free surface.Such hampering of the VEGF diffusion may cause the VEGF concentration in the biogel surface to be lower in regions with attached cells and indirect contact with the medium.It may delay and hinder the formation of a cellular network, which should be addressed in future work.
In the era when VEGF functions were unclear, theoretical works assumed that cellular networks developed purely due to the biogel deformation caused by the traction of endothelial cells without including bound VEGF [27][28][29].Such theoretical approaches obtained similar co-localized cells and biogel (Fig 2).We also considered gel deformation under the friction between the gel and the Petri dish so that our results could better explain the tilted biogel experiments [43].Comparing our results to the more recent mechanical-based model from van Oers et al. [37], we neglected the effects of stress gradients inside the cells.van Oers et al. considered that the intense substrate strain directly guides cell movement and deformation.They showed this mechanism can reproduce network formation from scattered cells and sprouting from cell spheroids.Since we studied cell morphology change signaled by soluble VEGF, we thus neglected substrate stress gradients inside the cell regions.The substrate deformation in the present model induced the bound VEGF gradients, which guided the cell directional migration.The soluble and bound VEGF effects on cell morphology and migration should be regarded as a feature for endothelial cells rather than for all the attached-type cell lines [58].The bound VEGF was consensusly regarded as dominating cell migration and network formation in the chemical-stimulation part.The oriented cue was sourced from the nonuniform distributed binding on biogel [35,42,59] by further incorporating the soluble VEGF in the model.

Conclusion
We have developed a comprehensive model encompassing mechanical regulation and chemical stimulation cues for forming capillary-like networks.The model covers the respective mechanisms of two VEGF forms on endothelial cells.The soluble VEGF in the culture medium induces morphological changes, and the bound VEGF on the biogel surface provides oriental signals, both playing roles in forming the capillary-like networks.The simulation successfully demonstrated that cells gathered from scattered individuals to cell cords and then from the cell cords to an overall network.The results were consistent with the experiments, showing that the cellular network was well-developed if the mechanical and chemical coupling worked simultaneously.The current model is valid for short-term static incubation.Furthermore, for the practical transplantation into living bodies to promote neovascularization around the surgical site [18,51,60], an extension of our model may cover the extracellular remodeling and endothelial cell morphology in the future.

Methods
The theoretical formulation depicts the in vitro capillary network development, considering the interactions among the endothelial cells, the biogel substrate, and the VEGF in soluble and bound forms.A two-dimensional hybrid model accounts for the endothelial cell adhesion and migration on the biogel surface, neglecting the cell penetration into the gel [43,61].As the schematic diagram shown in Fig 12 , The model considers the biogel deformation due to traction by the attached cells.The un-uniform distribution of bound VEGF due to the gel deformation offers a directional migration cue for the cells via haptotaxis.The soluble VEGF activates the morphology change so that the cells elongate in shape.As mentioned above, our mechanical mechanism is not direct between the cells and substrate.Still, it refers to the traction of biogel by cells, which provides a haptotaxis cue through the bound VEGF.The cellbased CPM predicts cellular individual behaviors, including migration, shape deformation, and morphology change [31,32].The continuum-based PDEs formulate the temporal and spatial distributions of biogel and VEGFs [33][34][35]42].

Cellular Potts model
Cell proliferation and apoptosis are negligible for the simulation of short-term incubation [5,6].Based on the two-dimensional assumption, the hollowing structure in cell morphology is also neglected in our model [52].In CPM, each cell occupies several connected Cartesian lattices, and a hypothetical energy of the cell system is defined using the following Hamiltonian [31,33] The variable x presents a specific lattice site and x' its neighboring lattices.We considered the first-and second-order neighborhood in the present modeling (Fig 13).The subindex τ x = 0, 1 in Eq (6) denotes whether the lattice x is cell-free or cell-occupied, respectively.The subindex σ x = n, a non-negative integer, represents the lattice precisely occupied by the n-th cell and σ x = 0 if the lattice is cell-free.The first term in Eq (6) describes the interface energy summation of the cell system, in which the interface energy J t x t x 0 between the two neighboring lattices takes different values depending on the neighboring states, and a Kroneckor delta formulation ð1 À d s x ;s x 0 Þ excludes counting the interface energy for the two neighboring lattices belonging to the same cell.The second-order neighbors contribute less energy with the weight 1= ffi ffi ffi 2 p than the first-order neighbors, considering their distance is ffi ffi ffi 2 p times longer.The second term in Eq (6) accounts for the cell deformation limitations, where a σ , p σ and l σ are the cell area, perimeter, and head-to-tail length, respectively, and λ a , λ p , λ l the corresponding restraint coefficients.The first and second restraints keep the initial cell area a 0 and maintain a minimal peripheral size to prevent the cell from breaking into pieces [31,42].The third restraint guides the endothelial cells to undergo cytoskeletal remodeling by changing their shape from round to elongated.The index ϕ σ = 0,1 labels whether the cell is still in its initial flattened state (ϕ σ = 0) or has changed its morphology (ϕ σ = 1).Once the cell has activated its shape remodeling, it is keen to increase its length l σ , presenting the differentiation for endothelial cells as irreversible [54].
The cell migration simulation followed the modified Metropolis algorithm [62].First, a lattice x in the computational domain and one of its neighbors x' were chosen, where the picking process was random and independent of the current lattice state.The state of the lattice x was assumed to extend to the lattice x' and thus caused the energy variation ΔH in the cell system.The Hamiltonian variation was further modified to include ΔH m for the haptotaxis effect by considering the bound VEGF gradient [24,25] The coefficient μ h represents the haptotaxis strength, and the b's are the bound VEGF concentrations at the two lattices.The term ð1 À d t x t x 0 Þ excludes haptotactic migration if both sites are occupied by cells, accounting for the phenomenon of cell contact inhibition [63].Whether or not the displacement happens was determined by the following probability [31] PðdisplacementÞ ¼

1
; The Potts model commonly takes the form exp(−ΔH m /T), including an additional temperature parameter T to quantify the magnitude fluctuation [31].We kept the temperature parameter constant without considering the fluctuation effect.A random number between 0 and 1 was generated and compared with the probability to decide whether the state change at x' truly happens.Eq (8) allows a possible move even under energy-raising conditions, representing the random walks for living cells.Finally, by repeating the above steps until the number of repetitions reaches the total lattice number, each lattice would have an equal opportunity for its state update, and the system would go through a complete Monte-Carlo step (MCS) equivalently [64].
The current model assumed all the cells were initially round-shaped.The possibility of the cell morphological change at each time step depended on the local concentration of soluble VEGF, which takes the following form.
where μ m is the activation coefficient, c σ the soluble VEGF concentration averaged over the identical cell-occupied lattices, and r σ the no-cell-contact interface ratio to its total perimeter.As reported in literature experiments, we included r σ keeping the cell in a cobblestone-like group to maintain its round, flattened shape [7,44,45].

Biogel PDE model
We neglected the secretion of extracellular matrix (ECM) by the endothelial cells, the limited enzymatic degradation rate of biogel, and the biogel proteolysis, considering the short-term incubation [7,65].Therefore, the time change rate of the local gel density ρ due to the biogel displacement is described by the conservation equation [27] @r @t ¼ where u is the biogel displacement vector.The biogel is considered a viscoelastic substrate.By ignoring the inertia force compared to the visco-elastic force, we adopted the following forcebalanced equation to evaluate the biogel displacement [29] r The first term σ cell represents the cell traction stress exerted on the surface, which is assumed proportional to the biogel density and formulated as where the parameter κ regulates the traction amplitude, the index τ x = 0, 1 coupled with the CPM indicates the stress existing only at the cell locations and I is the identity tensor.The second term σ vis in Eq (11) represents the viscous stress of the biogel.According to the Stokes law of deformation, the viscous stress is modeled as [27] Eq (13) assumes the viscous stress is proportional to the biogel strain and dilatation rates.In the equation ε ¼ ðru þ ðruÞ T Þ=2 is the gel strain, θ = r�u the gel dilatation, μ 1 and μ 2 the shear and bulk viscosities, respectively.The terms σ ela,linear and σ ela,long-range in Eq (11) are the local and long-range linear elastic stresses, respectively.In addition to the local stress, we assume that the elastic force includes that from the far-apart surroundings and is transmitted via the fiber components in the biogel.The local linear elastic stress takes the following form following Hooke's law [27] σ ela;linear ¼ where E is the Young's modulus and υ is the Poisson ratio of the biogel.The long-range elastic stress is modeled with the Laplacian terms for its ranged effects [29]: where β 1 and β 2 are the non-negative long-range elastic coefficients.The last term R ext in Eq (11) represents the friction force on the biogel by the Petri dish.As detailed in S1 Appendix, the attachment effect is expressed by the following form where h is the biogel thickness, which depends on the local gel density ρ, the initial thickness h 0 , and the initial density ρ 0

VEGFs PDE model
In general, the VEGF dynamics in cellular incubation systems contain the secretion from specific cell types, cell consumption via pinocytosis, and molecular degradation in the culture medium [23,42].For simplicity, we assumed the total amount of VEGF was a constant for simulating the short-term cultivation in this work.We considered two VEGF forms in the current model.The soluble form, the solute VEGF in the culture medium, can stimulate the cell morphological remodeling of elongation in cell shape.The bound form, the VEGF binding to biogel, can induce cell migration toward the haptotaxis cues [24,25,44].The two VEGF forms transform mutually via binding sites such as fibronectin and heparin sulfates in the biogel.Following the experimental designs in the literature [25,65], we presumed a finite amount of fibronectin in the biogel.The dynamic balance between the local concentrations of the soluble VEGF c and the bound VEGF b were formulated as where D c is the diffusion coefficient of soluble VEGF, k on and k off respectively represent the binding and unbinding rates of the VEGFs, and γ is the content ratio of VEGF binding sites to the biogel, wherein we assumed one fibronectin molecule provides a binding site [14].Finally, the term (γρ−bM FN /M VEGF ) evaluates the concentration of the free binding sites that remain conformed to the aforementioned binding constraint by including the molecular weight ratio of fibronectin M FN to VEGF M VEGF .

Numerical simulation
The hybrid model includes the CPM and PDEs.In-house coding with MATLAB was for CPM computation.COMSOL Multiphysics was adopted to solve the PDEs.The computational domain spanned 170x170 lattices.We set the squared lattice as 10 μm in width and MCS as 2 minutes so that the geometric center of each cell was monitored migrating with a maximum 30 μm/hr speed [66].Initially, a cell in a virtually round shape occupied 21 connected lattices (Fig 13).At the beginning of the simulation, 400 cells spread randomly in the lattice grid so that the number density for living cells was spatially comparable with the experimental reference [44].The grid size for the PDE models was identical to that of the CPM lattice.The biogel density and soluble VEGF concentration initially had a uniform distribution, whereas the initial concentration of bound VEGF was zero.Periodic conditions were specified at the four boundaries for the CPM and PDEs.These boundaries allow cells to migrate across and show up on the opposite side equivalently.
Table 1 presents the model parameters and their typical values.The endothelial cell line in the simulation referred to the human umbilical vein endothelial cells (HUVECs), and the primary biogel composition was the air-dried type I collagen.We determined the "This work" values in Table 1 by fitting the cell patterns with previous experiment data [44].The morphology coefficient μ m was first determined to have approximately one-third of the cells become elongated after 12 hours of incubation, which was the observed amount from the reference experiments [44].The cell-shape-related parameters λ a and λ p were then determined based on the round cell features of approximately 50 μm in diameter and 1960 μm 2 in area, and λ l was determined based on the elongated cell feature of 2.5 times the initial diameter [44].For the biogel parameters, the traction amplitude κ was determined so that the biogel experienced deformation comparable with previous simulation works [29,59], and the amount of fibronectin was chosen via the fibronectin ratio γ, which allowed approximately 10% of bound VEGF to transform from the soluble form [67]. Finally, the haptotaxis coefficient μ h was determined so that the overall cell pattern approached a capillary-like network after 18 hours of incubation, as shown in the experiments [7,43,44].

Images and statistical analysis
In this work, the number of cellular network junctions, segments, and the total length of the segments in the computational domain were adopted as indicators to quantify the cellular network integrity and complexity [14,68].A 100 μm threshold was applied to identify the network segments, and the network nodes intersected by at least three segments were judged as junctions.We further performed the cell shape analysis with the roundness index defined as: The dip unimodal test was performed for the cell migration trends [53].We adopted MATLAB R2020b and COMSOL Multiphysics 5.6 for numerical simulations, Image J 1.52a with the Angiogenesis Analyzer plugin for the network completeness check and the cell shape studies, and Microsoft Excel 2016 for Student's t-test statistical analysis.The regression analysis for the junction number and the dip tests were also performed using MATLAB R2020b.

Fig 1 .
Fig 1.Comparison between the literature experimental and current simulative cellular networks.(a) The experimental result was obtained at 18 hours of incubation of HUVECs [44].(b) The simulative cell distribution at 540 MCS (18 hours).Both diagrams show the cords connected by elongated cells and the cobblestone-like alliances composed of flattened round cells.https://doi.org/10.1371/journal.pcbi.1012281.g001

Fig 2 .
Fig 2. Evolution of the cell distribution, biogel density, and bound VEGF concentration with incubation time.The initial distributions were shown as 2 minutes.(a) The cell cords formed at 6 hours, connected sustainably during 6 to 18 hours and stabilized between 18 and 24 hours.(b) The biogel was locally concentrated by cell traction.(c) Bound VEGF on the biogel was concentrated following the biogel deformation.https://doi.org/10.1371/journal.pcbi.1012281.g002

Fig 3
shows the cell distributions, in which significant differences appear among the four simulation conditions.Benefiting from the periodic boundary settings for simulation, we could periodically arrange four identical computed cell images to become a 4-times figure, which could more effectively show the complete network lacunae surrounded by cellular cords than the original single figure.Unlike the typical case with both the cellular morphology and haptotaxis turned on (Fig 3A), the cellular cords composed of flattened, round cells could hardly connect well to develop into a complete network pattern if we turned off the cell morphology effect induced by soluble VEGF as shown by the two cases with the cellular morphology-off (Fig 3B and 3D).In the two cases with haptotaxis-off (Fig 3C and 3D), cell migration would only primarily follow random walks because of no spatial cues from the bound VEGF.Because the elongated cells could more easily touch each other than the round cells, a randomly constructed network formed in the case with morphology turned on (Fig 3C).If we further turned off the morphology effect, the cells were virtually in a spread distribution (Fig 3D).

Fig 3 .
Fig 3. Cell distribution pattern under morphology and haptotaxis effects turned on/off.The results show the complete network formed when morphology and haptotaxis effects were both on (case a).The un-connected cell cords occurred because of insufficient cord length in the two morphology-off cases (b and d).The randomly unsustainably connected cords happened in the two haptotaxis-off cases (c and d).https://doi.org/10.1371/journal.pcbi.1012281.g003

Fig 4 .
Fig 4. Junction number, segment number, and total segment length of the capillary-like network against incubation time.(a) Data from literature experiments [44].(b) Data from the current simulations under morphology and haptotaxis effects turned on/off.The simulations show networks were well-developed only when the morphology and haptotaxis mechanisms worked appropriately.The junction number, segment number, and total segment length in the two morphology-off cases were significantly lower than the two morphology-on cases.The junction number in the morphology-on/haptotaxis-off was considerably lower than in the morphology-on/haptotaxis-on.Each simulation case was repeated three times, and the data are presented as mean with standard deviation.The Student's t-test was performed towards the morphology-on and haptotaxis-on case: * = p < 0.05, ** = p < 0.01, and ns = non-significant.https://doi.org/10.1371/journal.pcbi.1012281.g004

Fig 5 .Fig 6 .Fig 7 .
Fig 5. Histogram of cell number versus (a) cell roundness and (b) cell aspect ratio at every 6 hours of incubation in the morphology-on/haptotaxis-on case.The cell roundness index, unity for an exact round shape, decreases when the cell elongates.The cell aspect ratio, unity for a round shape, increases when the cell elongates.The results show that the number of cells increased with time as they lost roundness, and the number of cells that kept aspect ratio of unity decreased with time.https://doi.org/10.1371/journal.pcbi.1012281.g005

Fig 8 .
Fig 8. Cell migration modality versus time under morphology and haptotaxis effects turned on/off.The classification was conducted on a two-hour basis among 400 cells.The results show that most cells performed random walks in the isotropic movement type for the four simulation cases.In the two haptotaxis-on cases (a and b), the onedirectional movement cells (red) were slightly more active in the first few hours.In the two morphology-on cases (a and c), the head-and-tail movement cells (blue) appear in the late hours.Each simulation was repeated three times for the same parameters.https://doi.org/10.1371/journal.pcbi.1012281.g008

Fig 9 .
Fig 9. Network junction number at 24 hours of incubations for various values of VEGF parameters.In each set, one parameter changes as denoted in the abscissa, and the other parameters keep their typical values as in Table1.The asterisk denotes the junction number at the typical parameter values of Table1.The network pattern for the 2x and 4x corresponding typical parameter values is shown on the right.The junction number in the network is denoted beside the network.A more complete network is of a more significant junction number, which approaches the asymptotic value as the parameter value increases, and there is virtually no open loop in sight when the junction number exceeds 170.

Fig 10 .
Fig 10.Temporal evolution of the junction number for different cell migration parameters.The results show that increasing these four parameter values raised the junction number in stable cellular networks.Increasing the initial VEGF concentration increased the junction number more significantly than the other parameters.In each figure, all the parameters are as the typical values in Table 1, except for the parameter in the box, which is adjusted by multiplying its typical value.https://doi.org/10.1371/journal.pcbi.1012281.g010

Fig 11 .
Fig 11.Relationships of the junction number at 24 hours of incubations with biogel parameters.The regression curves that fit the data trend depend on the data distribution within the simulation intervals.In each figure, one parameter changes as denoted in the abscissa, and the other parameters keep the typical values as in Table1.

Fig 12 .
Fig 12. Schematic diagram of (a) the cells-VEGF-biogel interplay and (b) the mechanisms included in the CPM and PDEs model.The cells exert traction force on substrates, causing the biogel deformation and non-uniform bound VEGF distribution.The bound VEGF gradient signal directed cell migration, resulting in a positive loop for cell aggregation, biogel deformation, and bound VEGF local accumulation.The two forms of VEGF transform dynamically.The soluble VEGF stimulates endothelial cells to initiate the elongated morphological change.https://doi.org/10.1371/journal.pcbi.1012281.g012

Fig 13 .
Fig 13.Initial cell shape and the n-th order neighborhoods in CPM.An individual cell occupies 21 connected lattices, initially approximating a round shape.The marked center lattice demonstrates four different types of neighboring lattices.The higher-order neighboring lattices are less likely to be selected as the extended site x' because they are far from the lattice x.We adopted the first and second-order neighborhoods in our simulation.The secondorder neighboring lattices are 1= ffi ffi ffi 2 p less likely to be selected than the first-order ones.https://doi.org/10.1371/journal.pcbi.1012281.g013

Table 1 .
The asterisk denotes the junction number at the typical parameter values of Table1.The network pattern for the 2x and 4x corresponding typical parameter values is shown on the right.The junction number in the network is denoted beside the network.A more complete network is of a more significant junction number, which approaches the asymptotic value as the parameter value increases, and there is virtually no open loop in sight when the junction number exceeds 170.https://doi.org/10.1371/journal.pcbi.1012281.g009